#### Analysis 
#### Jonne Kamphorst
#### April 2021 (origial file), February 2023 (updated figures for R&R). July 2023 (reproduction files finished)
#### R version: 4.3.0




#### Libraries and data ####
library(extrafont)
library(ggplot2)
library(estimatr)
library(margins)
library(ggeffects)
library(ggplot2)
library(viridis)
library(cregg)
library(dplyr)
library(ggpubr)
library(tidyverse)
library(cowplot)
library(mediation)

#font_import()
loadfonts(quiet = T, device = "pdf")
windowsFonts(Georgia = windowsFont("Georgia")) #load fonts for windows machines


## Set ggplot theme with the Georgia font
theme_set(theme_light(base_family = "Georgia", base_size = 12))

# ##Only keep the variables that I use in the analysis (JoP reproduction requirement), and anonymize the RIDs
#all_data <- readRDS("data/final_data_long_all_recoded.rds")
#all_data <- all_data %>% dplyr::select(rid, salience, partisan, partisan_altelect, issue, 
#                                       Partij, strategy_recoded, selected, position,
#                                       QPT2.1, QPT2.2, QPT2.3, QPT2.4, QPT3.1, QPT3.4, QPT3.5, task,
#                                       salience, QPT3.2)
#rids <- all_data %>% distinct(rid)
#random_numbers <- sample(1000:999999999999, size=nrow(rids), replace = FALSE)
#rids$random_numbers <- random_numbers
#all_data <- right_join(all_data, rids)
#all_data <- all_data %>% dplyr::select(-rid) %>% rename(rid = random_numbers)
#
#write_rds(file="data/final_data_long_all_recoded_reprod.rds", x=all_data)






#### Models for the long analysis #### 
all_data <- readRDS("data/final_data_long_all_recoded_reprod.rds")

## combine 4, 5 and 6 due to low number of people in those categories 
all_data$salience_recoded <- all_data$salience
all_data$salience_recoded[all_data$salience_recoded >4] <- 4
all_data$salience_recoded <- as.factor(all_data$salience_recoded)

#model as in PAP, w partisanship  
with_p <- (lm_robust(selected~strategy_recoded  * salience_recoded * partisan + issue + 
                       Partij * strategy_recoded, data=all_data, clusters = rid)) #add interaction term between attributes that were restricted

# model as in PAP yet without partisan, as in paper due to pre-treatment bias
without_p <- (lm_robust(selected~strategy_recoded * salience_recoded + issue + 
                          Partij * strategy_recoded, data=all_data, clusters = rid))

## Same models for the mechanisms
# model without partisan
without_p_pos <- (lm_robust(position~strategy_recoded * salience_recoded + issue + 
                              Partij * strategy_recoded, data=all_data, clusters = rid))










#### Plot the models: main models ####
## Marginal effects (differences in predicted probabilities) ##
saliences_recoded_levels <- as.factor(c(0, 1, 2, 3, 4)) 
partisan_recoded_levels <- as.factor(c(0, 1)) 

#run the margins (Takes a while, could be done faster if coded-up more efficiently, yet other packages don't give SEs)
#without_p_margins <- margins_summary(without_p, at=list(salience_recoded=saliences_recoded_levels))
#with_p_margins_1 <- margins_summary(with_p, at=list(salience_recoded=saliences_recoded_levels, partisan=1))
#with_p_margins_0 <- margins_summary(with_p, at=list(salience_recoded=saliences_recoded_levels, partisan=0))
with_p_margins <- full_join(with_p_margins_1, with_p_margins_0)







#Both margins without partisanship (Figure 1 left panel)
without_p_margins_plot_both <- ggplot(data=subset(without_p_margins, factor=="strategy_recodedDisagree" | 
                                                    factor =="strategy_recodedAmbiguity")) + 
  geom_pointrange(aes(y=AME, x=salience_recoded, ymin=lower, ymax=upper, color=factor),
                  position = position_dodge(width = 0.20), size=0.6) +
  labs(y ="Average Marginal Component Effect", x="Issue salience",
       title="Agreeing versus Disagreeing or Ambiguity",
       shape="Strategy",
       color="Strategy") + 
  geom_hline(yintercept = 0, linetype="dashed", color="black", size=1) + coord_cartesian(ylim = c(-0.3, 0.05)) +
  scale_colour_grey(end = 0.7, start=0.5, labels = c("Ambiguity", "Disagree")) +
  geom_line(aes(y=AME, x=salience_recoded, group=factor, color=factor), position = position_dodge(width = 0.20))
without_p_margins_plot_both
#ggsave("plots/fig1l_without_p_margins_plot_both.png", plot=without_p_margins_plot_both, dpi=600, width=3584, height=2616, units = c("px"))
#ggsave("plots/fig1l_without_p_margins_plot_both.tiff", plot=without_p_margins_plot_both, dpi=600, device="tiff", width=3584, height=2616, units = c("px"))



### With partisanship (Right panel Figure A3)
with_p_margins <- with_p_margins %>% rename(Strategy = factor,
                                            Partisanship = partisan)
with_p_margins$Strategy[with_p_margins$Strategy == "strategy_recodedAmbiguity"] <- "Ambiguity"
with_p_margins$Strategy[with_p_margins$Strategy == "strategy_recodedDisagree"] <- "Disagree"
with_p_margins$Partisanship[with_p_margins$Partisanship == 1] <- "Partisan"
with_p_margins$Partisanship[with_p_margins$Partisanship == 0] <- "Non-Partisan"
with_p_margins$salience_recoded <- as.numeric(with_p_margins$salience_recoded) #watch 0 to 5!



with_p_margins_plot_both <- ggplot(data=subset(with_p_margins, Strategy=="Ambiguity" | 
                     Strategy =="Disagree")) + 
  geom_pointrange(aes(y=AME, x=salience_recoded, ymin=lower, ymax=upper, color=Strategy, shape=Partisanship),
                  position = position_dodge(width = 0.20), size=0.6) +
  labs(y ="Average Marginal Component Effect", x="Issue salience",
       title="Agreeing versus Disagreeing or Ambiguity by Partisanship") +
       geom_hline(yintercept = 0, linetype="dashed", color="red", size=1) + coord_cartesian(ylim = c(-0.3, 0.05)) +
  scale_colour_grey(end = 0.6) +
  geom_line(aes(y=AME, x=salience_recoded, color=Strategy, linetype=Partisanship), position = position_dodge(width = 0.20))
with_p_margins_plot_both
#ggsave("plots/figa3r_with_p_margins_plot_both.png", plot=with_p_margins_plot_both, dpi=600, width = 6, height= 6)
#ggsave("plots/figa3r_with_p_margins_plot_both.tiff", plot=with_p_margins_plot_both, width = 6, height= 6, dpi=600, device="tiff")





rm(estimate, label, conf.low, conf.high, levels, plot_frame)







#### Predicted probabilities ####
without_p_predp <- ggeffect(without_p, terms=c("salience_recoded", "strategy_recoded")) #computes at the mean of a variable


# Predicted probability plots (right panel Figure 1)
without_p_predp_plot <- ggplot(data=without_p_predp, aes(x=x, y=predicted, ymin=conf.low, 
                                                         ymax=conf.high, color=group, shape=group, linetype=group)) +
  geom_pointrange(position = position_dodge(width = 0.10), size=0.7) + 
  labs(y ="Marginal Mean of Candidate Selection", x="Issue salience",
       title="Marginal Mean of Candidate Selection by Candidate Strategy",
       color="Strategy", shape="Strategy", group="Strategy", linetype="Strategy") + 
  geom_line(aes(y=predicted, x=x, group = group, color=group, linetype=group), size=0.8) +
  scale_color_grey(end=0.8) + coord_cartesian(ylim=c(0.345, 0.625))
without_p_predp_plot
#ggsave("plots/fig1r_without_p_predp_plot.png", plot=without_p_predp_plot, dpi=600, width=3584, height=2616, units = c("px"))
#ggsave("plots/fig1r_without_p_predp_plot.tiff", plot=without_p_predp_plot, dpi=600, device="tiff", width=3584, height=2616, units = c("px"))













#### Plot the models: mechanisms ####
## Marginal effects (differences in predicted probabilities) ##

#run the margins
#without_p_margins_pos <- margins_summary(without_p_pos, at=list(salience_recoded=saliences_recoded_levels))



### Figure 3. 
without_p_margins_plot_both_pos <- 
ggplot(data=subset(without_p_margins_pos, factor=="strategy_recodedDisagree" | 
                                                    factor =="strategy_recodedAmbiguity")) + 
  geom_pointrange(aes(y=AME, x=salience_recoded, ymin=lower, ymax=upper, color=factor),
                  position = position_dodge(width = 0.20), size=0.6) +
  labs(y ="Average Marginal Component Effect", x="Issue salience",
       title="Agreeing versus Disagreeing or Ambiguity",
       shape="Strategy",
       color="Strategy",
       subtitle="Outcome: `Does the candidate agree with you?'") + 
  geom_hline(yintercept = 0, linetype="dashed", color="black", size=1) +
  scale_colour_grey(end = 0.7, start=0.5, labels = c("Ambiguity", "Disagree")) +
  geom_line(aes(y=AME, x=salience_recoded, group=factor, color=factor), position = position_dodge(width = 0.20))
without_p_margins_plot_both_pos
ggsave("plots/fig3_without_p_margins_plot_both_pos.png", plot=without_p_margins_plot_both_pos, dpi=600, width=3584, height=2616, units = c("px"))
ggsave("plots/fig3_without_p_margins_plot_both_pos.tiff", plot=without_p_margins_plot_both_pos, dpi=600, device="tiff", width=3584, height=2616, units = c("px"))















#### Threshold value plots #####
without_p_predp <- ggeffect(without_p, terms=c("salience_recoded", "strategy_recoded"))



salience <- 0:4
noeff <- rep(0.5, 5)
amb <- without_p_predp$predicted[without_p_predp$group=="Ambiguity"]
agr <- without_p_predp$predicted[without_p_predp$group=="Agree"]
dis <- without_p_predp$predicted[without_p_predp$group=="Disagree"]

plot_frame <- data.frame(salience=salience, noeff=noeff,
                               amb=amb, agr=agr, dis=dis)
rm(salience, noeff, amb, agr, dis)


#calculate values
plot_frame$criticalvalue_estimate <- (plot_frame$amb - plot_frame$dis) / 
  (plot_frame$agr - plot_frame$dis) 


#plot the values: non-simulations
critval_plot <- ggplot(plot_frame, aes(x=salience, y=criticalvalue_estimate)) + geom_point(size=3) + geom_line(linewidth=0.8) +
  labs(x="Salience", y="Proportion of voters that need to agree with\na party for ambiguity to be inferior to agreeing") +
  theme(axis.title.y  = element_text(size = 10))    

critval_plot # Figure 2. 
ggsave("plots/fig2_critval_plot.png", plot=critval_plot, dpi=600, height=4.7, width=5) 
ggsave("plots/fig2_critval_plot.tiff", plot=critval_plot, dpi=600, height=4.7, width=5, device="tiff") 





#### Counter factual with the current electorate ####
# Select only the variables used (JoP requirement)
#noncj <- readRDS("data/final_data_noncj.rds")
#noncj <- noncj %>% dplyr::select(QPT1.1_4, QV6_abs, QPT1.1_2, QV3_abs, QPT1.1_5, QV2_abs, 
#                                 QPT1.1_3, QV4_abs, QPT1.1_1, QV1_abs, QPT1.1_6, QV5_abs, Q3.9, QV1 )
#                                
#write_rds(file="data/final_data_noncj_reprod.rds", x=noncj)
noncj <- readRDS("data/final_data_noncj_reprod.rds")





without_p_predp <- ggeffect(without_p, terms=c("salience_recoded", "strategy_recoded")) #computes at the mean of a variable
without_p_predp <- as.data.frame(without_p_predp) 




### For the effect of each issue and each strategy, what is the total effect in the electorate (based on the proportions of people who hold such an opinion in the electorate)
#for LGBTQ  #QPT1.1_4 is LGBTQ in qual and 6 in QVSR. 
agree_wstatem_lgbtq = without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 0] * prop.table(table(noncj$QPT1.1_4, noncj$QV6_abs))[1,1] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 1] * prop.table(table(noncj$QPT1.1_4, noncj$QV6_abs))[1,2] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 2] * prop.table(table(noncj$QPT1.1_4, noncj$QV6_abs))[1,3] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 3] * prop.table(table(noncj$QPT1.1_4, noncj$QV6_abs))[1,4] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 4] * (prop.table(table(noncj$QPT1.1_4, noncj$QV6_abs))[1,6] + 
                                                                                         prop.table(table(noncj$QPT1.1_4, noncj$QV6_abs))[1,7] + 
                                                                                         prop.table(table(noncj$QPT1.1_4, noncj$QV6_abs))[1,5]) +
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 0] * prop.table(table(noncj$QPT1.1_4, noncj$QV6_abs))[2,1] + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 1] *  prop.table(table(noncj$QPT1.1_4, noncj$QV6_abs))[2,2] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 2] *  prop.table(table(noncj$QPT1.1_4, noncj$QV6_abs))[2,3] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 3] *  prop.table(table(noncj$QPT1.1_4, noncj$QV6_abs))[2,4] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 4] *  (prop.table(table(noncj$QPT1.1_4, noncj$QV6_abs))[2,5] + 
                                                                                             prop.table(table(noncj$QPT1.1_4, noncj$QV6_abs))[2,6] + 
                                                                                             prop.table(table(noncj$QPT1.1_4, noncj$QV6_abs))[2,7])   

ambiguity_lgbtq = without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 0] * prop.table(table(noncj$QV6_abs))[1] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 1] * prop.table(table(noncj$QV6_abs))[2] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 2] * prop.table(table(noncj$QV6_abs))[3] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 3] * prop.table(table(noncj$QV6_abs))[4] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 4] * (prop.table(table(noncj$QV6_abs))[6] + 
                                                                                             prop.table(table(noncj$QV6_abs))[7] + 
                                                                                             prop.table(table(noncj$QV6_abs))[5]) 
  

disagree_wstatem_lgbtq = without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 0] * prop.table(table(noncj$QPT1.1_4, noncj$QV6_abs))[2,1] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 1] * prop.table(table(noncj$QPT1.1_4, noncj$QV6_abs))[2,2] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 2] * prop.table(table(noncj$QPT1.1_4, noncj$QV6_abs))[2,3] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 3] * prop.table(table(noncj$QPT1.1_4, noncj$QV6_abs))[2,4] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 4] *  (prop.table(table(noncj$QPT1.1_4, noncj$QV6_abs))[2,5] + 
                                                                                          prop.table(table(noncj$QPT1.1_4, noncj$QV6_abs))[2,6] + 
                                                                                          prop.table(table(noncj$QPT1.1_4, noncj$QV6_abs))[2,7]) + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 0] * prop.table(table(noncj$QPT1.1_4, noncj$QV6_abs))[1,1] + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 1] *  prop.table(table(noncj$QPT1.1_4, noncj$QV6_abs))[1,2] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 2] *  prop.table(table(noncj$QPT1.1_4, noncj$QV6_abs))[1,3] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 3] *  prop.table(table(noncj$QPT1.1_4, noncj$QV6_abs))[1,4] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 4] *  (prop.table(table(noncj$QPT1.1_4, noncj$QV6_abs))[1,5] + prop.table(table(noncj$QPT1.1_4, noncj$QV6_abs))[1,6] + prop.table(table(noncj$QPT1.1_4, noncj$QV6_abs))[1,7])


# for EU. EU is QV3 and  QPT1.1_2
agree_wstatem_eu = without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 0] * prop.table(table(noncj$QPT1.1_2, noncj$QV3_abs))[1,1] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 1] * prop.table(table(noncj$QPT1.1_2, noncj$QV3_abs))[1,2] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 2] * prop.table(table(noncj$QPT1.1_2, noncj$QV3_abs))[1,3] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 3] * prop.table(table(noncj$QPT1.1_2, noncj$QV3_abs))[1,4] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 4] * (prop.table(table(noncj$QPT1.1_2, noncj$QV3_abs))[1,5]+ prop.table(table(noncj$QPT1.1_2, noncj$QV3_abs))[1,6]) + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 0] * prop.table(table(noncj$QPT1.1_2, noncj$QV3_abs))[2,1] + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 1] *  prop.table(table(noncj$QPT1.1_2, noncj$QV3_abs))[2,2] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 2] *  prop.table(table(noncj$QPT1.1_2, noncj$QV3_abs))[2,3] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 3] *  prop.table(table(noncj$QPT1.1_2, noncj$QV3_abs))[2,4] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 4] *  (prop.table(table(noncj$QPT1.1_2, noncj$QV3_abs))[2,5] + prop.table(table(noncj$QPT1.1_2, noncj$QV3_abs))[2,6])

ambiguity_eu = without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 0] * prop.table(table(noncj$QV3_abs))[1] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 1] * prop.table(table(noncj$QV3_abs))[2] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 2] * prop.table(table(noncj$QV3_abs))[3] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 3] * prop.table(table(noncj$QV3_abs))[4] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 4] * (prop.table(table(noncj$QV3_abs))[5] + prop.table(table(noncj$QV3_abs))[6]) 

disagree_wstatem_eu = without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 0] * prop.table(table(noncj$QPT1.1_2, noncj$QV3_abs))[2,1] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 1] * prop.table(table(noncj$QPT1.1_2, noncj$QV3_abs))[2,2] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 2] * prop.table(table(noncj$QPT1.1_2, noncj$QV3_abs))[2,3] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 3] * prop.table(table(noncj$QPT1.1_2, noncj$QV3_abs))[2,4] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 4] * (prop.table(table(noncj$QPT1.1_2, noncj$QV3_abs))[2,5] + prop.table(table(noncj$QPT1.1_2, noncj$QV3_abs))[2,6]) + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 0] * prop.table(table(noncj$QPT1.1_2, noncj$QV3_abs))[1,1] + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 1] *  prop.table(table(noncj$QPT1.1_2, noncj$QV3_abs))[1,2] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 2] *  prop.table(table(noncj$QPT1.1_2, noncj$QV3_abs))[1,3] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 3] *  prop.table(table(noncj$QPT1.1_2, noncj$QV3_abs))[1,4] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 4] *  (prop.table(table(noncj$QPT1.1_2, noncj$QV3_abs))[1,5] + prop.table(table(noncj$QPT1.1_2, noncj$QV3_abs))[1,6])


#Minimum wage is QPT1.1_5 and QV 2
agree_wstatem_wage = without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 0] * prop.table(table(noncj$QPT1.1_5, noncj$QV2_abs))[1,1] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 1] * prop.table(table(noncj$QPT1.1_5, noncj$QV2_abs))[1,2] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 2] * prop.table(table(noncj$QPT1.1_5, noncj$QV2_abs))[1,3] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 3] * prop.table(table(noncj$QPT1.1_5, noncj$QV2_abs))[1,4] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 4] * (prop.table(table(noncj$QPT1.1_5, noncj$QV2_abs))[1,5] + prop.table(table(noncj$QPT1.1_5, noncj$QV2_abs))[1,6] + prop.table(table(noncj$QPT1.1_5, noncj$QV2_abs))[1,7]) + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 0] * prop.table(table(noncj$QPT1.1_5, noncj$QV2_abs))[2,1] + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 1] *  prop.table(table(noncj$QPT1.1_5, noncj$QV2_abs))[2,2] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 2] *  prop.table(table(noncj$QPT1.1_5, noncj$QV2_abs))[2,3] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 3] *  prop.table(table(noncj$QPT1.1_5, noncj$QV2_abs))[2,4] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 4] *  (prop.table(table(noncj$QPT1.1_5, noncj$QV2_abs))[2,5] + prop.table(table(noncj$QPT1.1_5, noncj$QV2_abs))[2,6] + prop.table(table(noncj$QPT1.1_5, noncj$QV2_abs))[2,7])

ambiguity_wage = without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 0] * prop.table(table(noncj$QV2_abs))[1] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 1] * prop.table(table(noncj$QV2_abs))[2] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 2] * prop.table(table(noncj$QV2_abs))[3] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 3] * prop.table(table(noncj$QV2_abs))[4] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 4] * (prop.table(table(noncj$QV2_abs))[5] + prop.table(table(noncj$QV2_abs))[6] + prop.table(table(noncj$QV2_abs))[7]) 

disagree_wstatem_wage = without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 0] * prop.table(table(noncj$QPT1.1_5, noncj$QV2_abs))[2,1] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 1] * prop.table(table(noncj$QPT1.1_5, noncj$QV2_abs))[2,2] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 2] * prop.table(table(noncj$QPT1.1_5, noncj$QV2_abs))[2,3] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 3] * prop.table(table(noncj$QPT1.1_5, noncj$QV2_abs))[2,4] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 4] * (prop.table(table(noncj$QPT1.1_5, noncj$QV2_abs))[2,5]  + prop.table(table(noncj$QPT1.1_5, noncj$QV2_abs))[2,6]  + prop.table(table(noncj$QPT1.1_5, noncj$QV2_abs))[2,7]) + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 0] * prop.table(table(noncj$QPT1.1_5, noncj$QV2_abs))[1,1] + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 1] *  prop.table(table(noncj$QPT1.1_5, noncj$QV2_abs))[1,2] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 2] *  prop.table(table(noncj$QPT1.1_5, noncj$QV2_abs))[1,3] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 3] *  prop.table(table(noncj$QPT1.1_5, noncj$QV2_abs))[1,4] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 4] *  (prop.table(table(noncj$QPT1.1_5, noncj$QV2_abs))[1,5] + prop.table(table(noncj$QPT1.1_5, noncj$QV2_abs))[1,6] + prop.table(table(noncj$QPT1.1_5, noncj$QV2_abs))[1,7])


#climate change QV_4 and QPT1.1_3
agree_wstatem_climate = without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 0] * prop.table(table(noncj$QPT1.1_3, noncj$QV4_abs))[1,1] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 1] * prop.table(table(noncj$QPT1.1_3, noncj$QV4_abs))[1,2] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 2] * prop.table(table(noncj$QPT1.1_3, noncj$QV4_abs))[1,3] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 3] * prop.table(table(noncj$QPT1.1_3, noncj$QV4_abs))[1,4] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 4] * (prop.table(table(noncj$QPT1.1_3, noncj$QV4_abs))[1,5] + prop.table(table(noncj$QPT1.1_3, noncj$QV4_abs))[1,6] + prop.table(table(noncj$QPT1.1_3, noncj$QV4_abs))[1,7]) + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 0] * prop.table(table(noncj$QPT1.1_3, noncj$QV4_abs))[2,1] + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 1] *  prop.table(table(noncj$QPT1.1_3, noncj$QV4_abs))[2,2] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 2] *  prop.table(table(noncj$QPT1.1_3, noncj$QV4_abs))[2,3] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 3] *  prop.table(table(noncj$QPT1.1_3, noncj$QV4_abs))[2,4] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 4] *  (prop.table(table(noncj$QPT1.1_3, noncj$QV4_abs))[2,5] + prop.table(table(noncj$QPT1.1_3, noncj$QV4_abs))[2,6] + prop.table(table(noncj$QPT1.1_3, noncj$QV4_abs))[2,7])

ambiguity_climate = without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 0] * prop.table(table(noncj$QV4_abs))[1] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 1] * prop.table(table(noncj$QV4_abs))[2] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 2] * prop.table(table(noncj$QV4_abs))[3] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 3] * prop.table(table(noncj$QV4_abs))[4] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 4] * (prop.table(table(noncj$QV4_abs))[5] + prop.table(table(noncj$QV4_abs))[6] + prop.table(table(noncj$QV4_abs))[7]) 

disagree_wstatem_climate = without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 0] * prop.table(table(noncj$QPT1.1_3, noncj$QV4_abs))[2,1] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 1] * prop.table(table(noncj$QPT1.1_3, noncj$QV4_abs))[2,2] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 2] * prop.table(table(noncj$QPT1.1_3, noncj$QV4_abs))[2,3] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 3] * prop.table(table(noncj$QPT1.1_3, noncj$QV4_abs))[2,4] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 4] * (prop.table(table(noncj$QPT1.1_3, noncj$QV4_abs))[2,5] + prop.table(table(noncj$QPT1.1_3, noncj$QV4_abs))[2,6] + prop.table(table(noncj$QPT1.1_3, noncj$QV4_abs))[2,7]) + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 0] * prop.table(table(noncj$QPT1.1_3, noncj$QV4_abs))[1,1] + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 1] *  prop.table(table(noncj$QPT1.1_3, noncj$QV4_abs))[1,2] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 2] *  prop.table(table(noncj$QPT1.1_3, noncj$QV4_abs))[1,3] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 3] *  prop.table(table(noncj$QPT1.1_3, noncj$QV4_abs))[1,4] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 4] *  (prop.table(table(noncj$QPT1.1_3, noncj$QV4_abs))[1,5] + prop.table(table(noncj$QPT1.1_3, noncj$QV4_abs))[1,6] + prop.table(table(noncj$QPT1.1_3, noncj$QV4_abs))[1,7])


#immigration QV_1 is and QPT1.1_1
agree_wstatem_immi = without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 0] * prop.table(table(noncj$QPT1.1_1, noncj$QV1_abs))[1,1] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 1] * prop.table(table(noncj$QPT1.1_1, noncj$QV1_abs))[1,2] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 2] * prop.table(table(noncj$QPT1.1_1, noncj$QV1_abs))[1,3] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 3] * prop.table(table(noncj$QPT1.1_1, noncj$QV1_abs))[1,4] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 4] * (prop.table(table(noncj$QPT1.1_1, noncj$QV1_abs))[1,5] + prop.table(table(noncj$QPT1.1_1, noncj$QV1_abs))[1,6] + prop.table(table(noncj$QPT1.1_1, noncj$QV1_abs))[1,7]) + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 0] * prop.table(table(noncj$QPT1.1_1, noncj$QV1_abs))[2,1] + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 1] *  prop.table(table(noncj$QPT1.1_1, noncj$QV1_abs))[2,2] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 2] *  prop.table(table(noncj$QPT1.1_1, noncj$QV1_abs))[2,3] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 3] *  prop.table(table(noncj$QPT1.1_1, noncj$QV1_abs))[2,4] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 4] *  (prop.table(table(noncj$QPT1.1_1, noncj$QV1_abs))[2,5] + prop.table(table(noncj$QPT1.1_1, noncj$QV1_abs))[2,6] + prop.table(table(noncj$QPT1.1_1, noncj$QV1_abs))[2,7])

ambiguity_immi = without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 0] * prop.table(table(noncj$QV1_abs))[1] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 1] * prop.table(table(noncj$QV1_abs))[2] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 2] * prop.table(table(noncj$QV1_abs))[3] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 3] * prop.table(table(noncj$QV1_abs))[4] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 4] * (prop.table(table(noncj$QV1_abs))[5] + prop.table(table(noncj$QV1_abs))[6] + prop.table(table(noncj$QV1_abs))[7]) 

disagree_wstatem_immi = without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 0] * prop.table(table(noncj$QPT1.1_1, noncj$QV1_abs))[2,1] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 1] * prop.table(table(noncj$QPT1.1_1, noncj$QV1_abs))[2,2] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 2] * prop.table(table(noncj$QPT1.1_1, noncj$QV1_abs))[2,3] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 3] * prop.table(table(noncj$QPT1.1_1, noncj$QV1_abs))[2,4] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 4] * (prop.table(table(noncj$QPT1.1_1, noncj$QV1_abs))[2,5]  + prop.table(table(noncj$QPT1.1_1, noncj$QV1_abs))[2,6] + prop.table(table(noncj$QPT1.1_1, noncj$QV1_abs))[2,7]) + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 0] * prop.table(table(noncj$QPT1.1_1, noncj$QV1_abs))[1,1] + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 1] *  prop.table(table(noncj$QPT1.1_1, noncj$QV1_abs))[1,2] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 2] *  prop.table(table(noncj$QPT1.1_1, noncj$QV1_abs))[1,3] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 3] *  prop.table(table(noncj$QPT1.1_1, noncj$QV1_abs))[1,4] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 4] *  (prop.table(table(noncj$QPT1.1_1, noncj$QV1_abs))[1,5]  + prop.table(table(noncj$QPT1.1_1, noncj$QV1_abs))[1,6] + prop.table(table(noncj$QPT1.1_1, noncj$QV1_abs))[1,7])



#farms QV_5 is farms and QPT1.1_6
agree_wstatem_farms = without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 0] * prop.table(table(noncj$QPT1.1_6, noncj$QV5_abs))[1,1] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 1] * prop.table(table(noncj$QPT1.1_6, noncj$QV5_abs))[1,2] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 2] * prop.table(table(noncj$QPT1.1_6, noncj$QV5_abs))[1,3] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 3] * prop.table(table(noncj$QPT1.1_6, noncj$QV5_abs))[1,4] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 4] * (prop.table(table(noncj$QPT1.1_6, noncj$QV5_abs))[1,5] + prop.table(table(noncj$QPT1.1_6, noncj$QV5_abs))[1,6] + prop.table(table(noncj$QPT1.1_6, noncj$QV5_abs))[1,7]) + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 0] * prop.table(table(noncj$QPT1.1_6, noncj$QV5_abs))[2,1] + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 1] *  prop.table(table(noncj$QPT1.1_6, noncj$QV5_abs))[2,2] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 2] *  prop.table(table(noncj$QPT1.1_6, noncj$QV5_abs))[2,3] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 3] *  prop.table(table(noncj$QPT1.1_6, noncj$QV5_abs))[2,4] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 4] *  (prop.table(table(noncj$QPT1.1_6, noncj$QV5_abs))[2,5]  + prop.table(table(noncj$QPT1.1_6, noncj$QV5_abs))[2,6] + prop.table(table(noncj$QPT1.1_6, noncj$QV5_abs))[2,7])

ambiguity_farms = without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 0] * prop.table(table(noncj$QV5_abs))[1] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 1] * prop.table(table(noncj$QV5_abs))[2] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 2] * prop.table(table(noncj$QV5_abs))[3] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 3] * prop.table(table(noncj$QV5_abs))[4] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 4] * (prop.table(table(noncj$QV5_abs))[5] + prop.table(table(noncj$QV5_abs))[6] + prop.table(table(noncj$QV5_abs))[7]) 

disagree_wstatem_farms = without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 0] * prop.table(table(noncj$QPT1.1_6, noncj$QV5_abs))[2,1] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 1] * prop.table(table(noncj$QPT1.1_6, noncj$QV5_abs))[2,2] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 2] * prop.table(table(noncj$QPT1.1_6, noncj$QV5_abs))[2,3] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 3] * prop.table(table(noncj$QPT1.1_6, noncj$QV5_abs))[2,4] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 4] * (prop.table(table(noncj$QPT1.1_6, noncj$QV5_abs))[2,5] + prop.table(table(noncj$QPT1.1_6, noncj$QV5_abs))[2,6] + prop.table(table(noncj$QPT1.1_6, noncj$QV5_abs))[2,7]) + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 0] * prop.table(table(noncj$QPT1.1_6, noncj$QV5_abs))[1,1] + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 1] *  prop.table(table(noncj$QPT1.1_6, noncj$QV5_abs))[1,2] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 2] *  prop.table(table(noncj$QPT1.1_6, noncj$QV5_abs))[1,3] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 3] *  prop.table(table(noncj$QPT1.1_6, noncj$QV5_abs))[1,4] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 4] *  (prop.table(table(noncj$QPT1.1_6, noncj$QV5_abs))[1,5]  + prop.table(table(noncj$QPT1.1_6, noncj$QV5_abs))[1,6] + prop.table(table(noncj$QPT1.1_6, noncj$QV5_abs))[1,7])



estimates <- c(agree_wstatem_climate, agree_wstatem_eu, agree_wstatem_farms, agree_wstatem_immi, agree_wstatem_lgbtq, agree_wstatem_wage,
            ambiguity_climate, ambiguity_eu, ambiguity_farms, ambiguity_immi, ambiguity_lgbtq, ambiguity_wage,
            disagree_wstatem_climate, disagree_wstatem_eu, disagree_wstatem_farms, disagree_wstatem_immi, disagree_wstatem_lgbtq, disagree_wstatem_wage)
labels <- c("Climate Pro", "EU Pro", "Farms Pro", "Immigration Pro","LGBTQ+ Pro", "Min. Wage Pro", 
            "Climate Ambiguity", "EU Ambiguity", "Farms Ambiguity", "Immigration Ambiguity", "LGBTQ+ Ambiguity", "Min. Wage Ambiguity",
            "Climate Anti", "EU Anti", "Farms Anti", "Immigration Anti", "LGBTQ+ Anti", "Min. Wage Anti")
Estimate <- c("Climate", "EU", "Farms", "Immigration", "LGBTQ+", "Min. Wage",
              "Climate", "EU", "Farms", "Immigration", "LGBTQ+", "Min. Wage",
              "Climate", "EU", "Farms", "Immigration", "LGBTQ+", "Min. Wage")

ASE <- mean(without_p_predp$std.error[without_p_predp$group=="Agree"])
AMBSE <- mean(without_p_predp$std.error[without_p_predp$group=="Ambiguity"])
DISSE <- mean(without_p_predp$std.error[without_p_predp$group=="Disagree"])
SE <- c(rep(ASE, 6), rep(AMBSE, 6), rep(DISSE, 6))

plot_frame_simulation <- data.frame(estimate = estimates, label = labels, Estimate = Estimate, SE= SE)
rm(estimates, labels, Estimate, ASE, AMBSE, DISSE, SE)


levels <- c("LGBTQ+ Anti", "LGBTQ+ Ambiguity", "LGBTQ+ Pro",
            "Immigration Anti", "Immigration Ambiguity", "Immigration Pro",
            "Min. Wage Anti", "Min. Wage Ambiguity", "Min. Wage Pro",
            "Farms Anti", "Farms Ambiguity", "Farms Pro",
            "Climate Anti", "Climate Ambiguity", "Climate Pro",
            "EU Anti", "EU Ambiguity", "EU Pro",
            "PvdA", "PVV", "VVD", "SP", "GL", "CDA")

plot_frame_simulation$label <- factor(plot_frame_simulation$label, levels = levels)

plot_frame_simulation$conf.high <- plot_frame_simulation$estimate + 1.96 * plot_frame_simulation$SE
plot_frame_simulation$conf.low <- plot_frame_simulation$estimate - 1.96 * plot_frame_simulation$SE


### Figure 5 left panel
simulation_counterfactuals <- ggplot(data=plot_frame_simulation, aes(x=estimate, y=label, 
                                                                     xmin=conf.low, xmax=conf.high)) +
  geom_pointrange(size=0.6) + 
  ylab("") + xlab("Mean probability of selection") + scale_colour_grey(end=0.8) +
  theme(legend.position = "none", strip.text = element_text(colour = 'black'),
        strip.background =element_rect(fill="lightgrey")) + labs(subtitle = "The effect of being Pro, Anti, or Ambiguous\non the policy statement") +
  facet_wrap(~Estimate, strip.position = "right", ncol=1, scales = "free_y")

simulation_counterfactuals
#ggsave("plots/fig5l_simulation_counterfactuals.png", plot=simulation_counterfactuals, dpi=600, width=5.5, height=7)
#ggsave("plots/fig5l_simulation_counterfactuals.tiff", plot=simulation_counterfactuals, dpi=600, device="tiff", width=5.5, height=7)








##### Add a simulation plot with the current left electorate for the PvdA ##### 517
# I.e. people who think favorable about the PvdA. Current best strategy on immigration, what if positions change? 

### Current PvdA electorate on immigration: show that now agreeing is better than disagreeing. Ambiguity best
# Add if someone is a pvda sympathized
noncj$therm_pvda <- noncj$Q3.9 
noncj$therm_pvda[noncj$therm_pvda == "Zeer sympathiek: 10"] <- "10"
noncj$therm_pvda[noncj$therm_pvda == "Zeer onsympathiek: 0"] <- "0"
noncj$therm_pvda[noncj$therm_pvda == "Weet ik niet"] <- NA
noncj$therm_pvda <- as.numeric(noncj$therm_pvda)

#only keep pvda sympathizers
noncj_pvda <- noncj %>% filter(therm_pvda > 6)


#immigration QV_1 is and QPT1.1_1
agree_wstatem_immi = without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 0] * prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs))[1,1] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 1] * prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs))[1,2] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 2] * prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs))[1,3] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 3] * prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs))[1,4] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 4] * (prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs))[1,6] + prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs))[1,5] + prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs))[1,7]) + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 0] * prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs))[2,1] + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 1] *  prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs))[2,2] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 2] *  prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs))[2,3] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 3] *  prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs))[2,4] + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 4] *  (prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs))[2,5] + prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs))[2,6]+ prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs))[2,7])

ambiguity_immi = without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 0] * prop.table(table(noncj_pvda$QV1_abs))[1] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 1] * prop.table(table(noncj_pvda$QV1_abs))[2] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 2] * prop.table(table(noncj_pvda$QV1_abs))[3] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 3] * prop.table(table(noncj_pvda$QV1_abs))[4] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 4] * (prop.table(table(noncj_pvda$QV1_abs))[5] + prop.table(table(noncj_pvda$QV1_abs))[6] + prop.table(table(noncj_pvda$QV1_abs))[7]) 

disagree_wstatem_immi = without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 0] * prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs))[2,1] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 1] * prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs))[2,2] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 2] * prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs))[2,3] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 3] * prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs))[2,4] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 4] * (prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs))[2,5] + prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs))[2,6] + prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs))[2,7]) + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 0] * prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs))[1,1] + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 1] *  prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs))[1,2] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 2] *  prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs))[1,3] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 3] *  prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs))[1,4] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 4] *  (prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs))[1,5] + prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs))[1,6] + prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs))[1,7])




### Now show what happens if progressives become more salient (thus positions stay the same)
## add a new salience variable. Add +2 for everyone who is eens, set max to 5. 

noncj_pvda$QV1_alt <- noncj_pvda$QV1
noncj_pvda$QV1_alt[noncj_pvda$QPT1.1_1 == "Eens"] <- noncj_pvda$QV1_alt[noncj_pvda$QPT1.1_1 == "Eens"] + 2
noncj_pvda$QV1_alt[noncj_pvda$QV1_alt>6] <- 6

noncj_pvda$QV1_abs_alt <- noncj_pvda$QV1_abs
noncj_pvda$QV1_abs_alt[noncj_pvda$QPT1.1_1 == "Eens"] <- noncj_pvda$QV1_abs_alt[noncj_pvda$QPT1.1_1 == "Eens"] + 2
noncj_pvda$QV1_abs_alt[noncj_pvda$QV1_abs_alt>6] <- 5

agree_wstatem_immi_alt = without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 0] * prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs_alt))[1,1] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 1] * prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs_alt))[1,2] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 2] * prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs_alt))[1,3] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 3] * prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs_alt))[1,4] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 4] * (prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs_alt))[1,5] + prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs_alt))[1,6] + prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs_alt))[1,7]) + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 0] * prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs_alt))[2,1] + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 1] *  prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs_alt))[2,2] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 2] *  prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs_alt))[2,3] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 3] *  prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs_alt))[2,4] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 4] *  (prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs_alt))[2,5] + prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs_alt))[2,6] + prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs_alt))[2,7])

ambiguity_immi_alt = without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 0] * prop.table(table(noncj_pvda$QV1_abs_alt))[1] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 1] * prop.table(table(noncj_pvda$QV1_abs_alt))[2] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 2] * prop.table(table(noncj_pvda$QV1_abs_alt))[3] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 3] * prop.table(table(noncj_pvda$QV1_abs_alt))[4] + 
  without_p_predp$predicted[without_p_predp$group=="Ambiguity" & without_p_predp$x== 4] * (prop.table(table(noncj_pvda$QV1_abs_alt))[5]  + prop.table(table(noncj_pvda$QV1_abs_alt))[6] + prop.table(table(noncj_pvda$QV1_abs_alt))[7] ) 

disagree_wstatem_immi_alt = without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 0] * prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs_alt))[2,1] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 1] * prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs_alt))[2,2] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 2] * prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs_alt))[2,3] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 3] * prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs_alt))[2,4] + 
  without_p_predp$predicted[without_p_predp$group=="Agree" & without_p_predp$x== 4] * (prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs_alt))[2,5] + prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs_alt))[2,6] + prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs_alt))[2,7]) + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 0] * prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs_alt))[1,1] + 
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 1] *  prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs_alt))[1,2] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 2] *  prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs_alt))[1,3] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 3] *  prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs_alt))[1,4] +  
  without_p_predp$predicted[without_p_predp$group=="Disagree" & without_p_predp$x== 4] *  (prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs_alt))[1,5] + prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs_alt))[1,6] + prop.table(table(noncj_pvda$QPT1.1_1, noncj_pvda$QV1_abs_alt))[1,7])








estimates <- c(agree_wstatem_immi, ambiguity_immi, disagree_wstatem_immi, agree_wstatem_immi_alt, ambiguity_immi_alt, disagree_wstatem_immi_alt)
labels <- c("Pro", "Ambiguity", "Anti", 
            "Pro","Ambiguity", "Anti")
Estimate <- c("Actual distribution", "Actual distribution", "Actual distribution",
              "Hypothetical distribution", "Hypothetical distribution", "Hypothetical distribution")

ASE <- mean(without_p_predp$std.error[without_p_predp$group=="Agree"])
AMBSE <- mean(without_p_predp$std.error[without_p_predp$group=="Ambiguity"])
DISSE <- mean(without_p_predp$std.error[without_p_predp$group=="Disagree"])
SE <- c(ASE, AMBSE, DISSE, ASE, AMBSE, DISSE)

plot_frame_simulation <- data.frame(estimate = estimates, label = labels, Estimate = Estimate, SE= SE)
rm(estimates, labels, Estimate, ASE, AMBSE, DISSE, SE)

plot_frame_simulation$conf.high <- plot_frame_simulation$estimate + 1.96 * plot_frame_simulation$SE
plot_frame_simulation$conf.low <- plot_frame_simulation$estimate - 1.96 * plot_frame_simulation$SE
plot_frame_simulation$label <- factor(plot_frame_simulation$label, levels=c("Pro", "Ambiguity", "Anti"))

# Plot the strategies
strategies <- ggplot(data=plot_frame_simulation, aes(y=estimate, x=label, ymin=conf.low, ymax=conf.high)) +
  geom_pointrange(size=0.6) + 
  facet_grid(cols=vars(Estimate), switch="x") +
  coord_flip() +  # flip coordinates (puts labels on y axis)
  xlab("") + ylab("Mean Pr(selection)") +
  theme(legend.position = "none",
        strip.text = element_text(colour = 'black'),
        strip.background =element_rect(fill="lightgrey")) + ggtitle("")



left <- ggplot(data=noncj_pvda, aes(x=QV1)) + 
  geom_density(size=1, bw=0.7) + theme_minimal(base_family = "Georgia") +
  theme(axis.text.y=element_blank(), 
        legend.position = "none") + 
  scale_x_continuous(breaks=seq(-6,6,1)) + xlim(c(-6,6)) + labs(y="Density", x="Salience (votes)") +
  scale_fill_grey(end=0.8) + ggtitle("", subtitle="")

left1 <- ggplot(data=noncj_pvda, aes(x=QV1_alt)) + 
  geom_density(size=1, bw=0.7) + theme_minimal(base_family = "Georgia") +
  theme(axis.text.y=element_blank(), 
        legend.position = "none") + 
  scale_x_continuous(breaks=seq(-6,6,1)) + xlim(c(-6,6)) + labs(y="", x="") +
  scale_fill_grey(end=0.8) + ggtitle("", subtitle="")


#### Add to one plot (Figure 6 - final placement of plots done post R)
library("gridExtra")
library("cowplot")
lay <- rbind(c(1,2),
             c(5,5),
             c(5,5))
PvdA_plots <- grid.arrange(left, left1,
                                  strategies, layout_matrix = lay)

PvdA_plots # Figure 6 
ggsave("plots/fig6_PvdA_plots.png", plot=PvdA_plots, dpi=600, width = 7, height = 5.2)
ggsave("plots/fig6_PvdA_plots.tiff", plot=PvdA_plots, dpi=600, width = 7, height = 5.2, device="tiff")











#### Add methods for direct vs indirect effects ####
# Fit the models
fit.mediator <- (lm(position ~ strategy_recoded * salience_recoded + issue, data=all_data, na.action=na.omit)) #number of observations must match between mediator and outcome models

fit.dv <- (lm(selected ~ strategy_recoded * salience_recoded + position + issue, 
              data=all_data, na.action=na.omit))

cluster <- all_data %>% drop_na(selected, strategy_recoded, salience_recoded, position, issue) %>%
  dplyr::select(rid) #cluster variable needs to be its own separate variable and not one in the dataset


## loop that gets the results for different levels of the moderating variable and puts these in a dataframe

#create empty objects
ACME <- rep(NA, length(0:4))
ADE <- rep(NA, length(0:4))
ACME_lo <- rep(NA, length(0:4))
ACME_hi <- rep(NA, length(0:4))
ADE_lo <- rep(NA, length(0:4))
ADE_hi <- rep(NA, length(0:4))
salience <- rep(NA, length(0:4))


#loop through the different potential levels of the salience variable
for (i in 0:4) {
  results <- mediate(fit.mediator, fit.dv, treat="strategy_recoded", mediator= "position", 
                   covariates = list(salience_recoded = i),
                   treat.value = "Ambiguity", control.value = "Agree", #same as in the other models
                   sims = 500,
                   cluster = cluster) 

  summary(results)
  ACME[1+i] <- results$d0
  ADE[1+i] <- results$z0
  ACME_lo[1+i] <- results$d0.ci[1]
  ACME_hi[1+i] <- results$d0.ci[2]
  ADE_lo[1+i] <- results$z0.ci[1]
  ADE_hi[1+i] <- results$z0.ci[2]
  salience[1+i] <- i
}

mediation_results <- data.frame(ACME=ACME, 
                         ADE=ADE,
                      ACME_lo=ACME_lo,
                      ACME_hi=ACME_hi,
                      ADE_lo=ADE_lo,
                      ADE_hi=ADE_hi,
                      salience=salience)

#change to dataframe for plotting
effects <- c(mediation_results$ACME, mediation_results$ADE)
ci.lo <- c(mediation_results$ACME_lo, mediation_results$ADE_lo)
ci.hi <- c(mediation_results$ACME_hi, mediation_results$ADE_hi)
estimate <- c(rep("ACME", 5), rep("ADE", 5))
salience <- c(0,1,2,3,4,0,1,2,3,4)

mediation_results_plot <- data.frame(
  effects=effects,
  ci.lo=ci.lo,
  ci.hi=ci.hi,
  estimate=estimate,
  salience=salience
)

#plot
mediation_plot <- ggplot(mediation_results_plot, aes(x=salience, y=effects, color=estimate, ymin=ci.lo, ymax=ci.hi)) +
  geom_pointrange(size=0.6, position = position_dodge(width = 0.2)) + 
  geom_hline(yintercept = 0, linetype="dashed", color="black", size=1) + 
  scale_color_grey(end=0.65) +
  labs(y="Average Marginaal Component Effect", x= "Issue salience", color="Estimate",
       title="")

mediation_plot


# model without the level of the covariate defined, this gives us the average proportion mediated as reported in the paper
# Prop. Mediated  0.24578      0.17802         0.32  <2e-16 ***
summary(mediate(fit.mediator, fit.dv, treat="strategy_recoded", mediator= "position", 
                treat.value = "Ambiguity", control.value = "Agree", #same as in the other models
                sims = 500,
                cluster = cluster)) 


### Figure 4.
ggsave("plots/fig4_mediation_plot.png", plot= mediation_plot, dpi=600, width=3584, height=2616, units = c("px"))
ggsave("plots/fig4_mediation_plot.tiff", plot= mediation_plot, dpi=600, device="tiff", width=3584, height=2616, units = c("px"))










#### Checking whether different salience groups are balanced ####
## Run a multinomial logit model in r and stata. Predictors are continues
##    Outcome is a factorized variable for the different salience levels. 
## Variables:
# age QPT2.1
# gender QPT2.2
# education QPT2.3
# income QPT2.4
# interest in pol QPT3.1
# competence to participate QPT3.2
# left right QPT3.4
# risk attitudes QPT3.5

all_data_forstata <- all_data %>% filter(task==1) 


all_data_forstata$age <- NA
all_data_forstata$age[all_data_forstata$QPT2.1 == "16-24"] <- 0
all_data_forstata$age[all_data_forstata$QPT2.1 == "25-34"] <- 1
all_data_forstata$age[all_data_forstata$QPT2.1 == "35-44"] <- 2
all_data_forstata$age[all_data_forstata$QPT2.1 == "45-54"] <- 3
all_data_forstata$age[all_data_forstata$QPT2.1 == "55-64"] <- 4
all_data_forstata$age[all_data_forstata$QPT2.1 == "65+"] <- 5

all_data_forstata$gender <- all_data_forstata$QPT2.2

all_data_forstata$education <- NA
all_data_forstata$education[all_data_forstata$QPT2.3 == "Basisschool"] <- 0
all_data_forstata$education[all_data_forstata$QPT2.3 == "LBO, VBO, LEAO, LTS ambachtsschool, huishoudschool, LHNO, VMBO (niveau 1-3), KADER."] <- 1
all_data_forstata$education[all_data_forstata$QPT2.3 == "MBO niveau 1 afgemaakt" ] <- 2
all_data_forstata$education[all_data_forstata$QPT2.3 == "MULO, ULO, MAVO, VMBO (niveau 4; theoretische leerweg), HAVO jaar 3-4; VWO jaar 3-5" ] <- 3
all_data_forstata$education[all_data_forstata$QPT2.3 == "MBO niveau 2 en 3 afgemaakt" ] <- 4
all_data_forstata$education[all_data_forstata$QPT2.3 == "KMBO, leerlingwezen, MBO, MEAO, MTS afgemaakt" ] <- 5
all_data_forstata$education[all_data_forstata$QPT2.3 == "MBO niveau 4 afgemaakt" ] <- 6
all_data_forstata$education[all_data_forstata$QPT2.3 == "HAVO, MMS, MSVM afgemaakt"] <- 7
all_data_forstata$education[all_data_forstata$QPT2.3 == "VWO, HBS, atheneum, gymnasium afgemaakt" ] <- 8
all_data_forstata$education[all_data_forstata$QPT2.3 == "Bachelor HBO, kweekschool, PABO, conservatorium, MO-akten afgemaakt"] <- 9
all_data_forstata$education[all_data_forstata$QPT2.3 == "Bachelor universiteit afgemaakt"] <- 10
all_data_forstata$education[all_data_forstata$QPT2.3 == "WO/universiteit: Master's degree, tweede fase opleidingen, ingenieur, meester, doctorandus"] <- 11
all_data_forstata$education[all_data_forstata$QPT2.3 == "Doctoraat / gepromoveerd"] <- 12
table(all_data_forstata$education)

all_data_forstata$income <- NA
all_data_forstata$income[all_data_forstata$QPT2.4 == "Weet ik niet of wil ik niet zeggen"] <- NA
all_data_forstata$income[all_data_forstata$QPT2.4 == "Minder dan â‚¬1.100"] <- 0
all_data_forstata$income[all_data_forstata$QPT2.4 == "â‚¬1.100 tot â‚¬1.400"] <- 1
all_data_forstata$income[all_data_forstata$QPT2.4 == "â‚¬1.400 tot â‚¬1.700"] <- 2
all_data_forstata$income[all_data_forstata$QPT2.4 == "â‚¬1.700 tot â‚¬2.000"] <- 3
all_data_forstata$income[all_data_forstata$QPT2.4 == "â‚¬2.000 tot â‚¬2.400"] <- 4
all_data_forstata$income[all_data_forstata$QPT2.4 == "â‚¬2.400 tot â‚¬2.790"] <- 5
all_data_forstata$income[all_data_forstata$QPT2.4 == "â‚¬2.800 tot â‚¬3.300"] <- 6
all_data_forstata$income[all_data_forstata$QPT2.4 == "â‚¬3.300 tot â‚¬3.900"] <- 7
all_data_forstata$income[all_data_forstata$QPT2.4 == "â‚¬3.900 tot â‚¬4.800"] <- 8
all_data_forstata$income[all_data_forstata$QPT2.4 == "â‚¬4.800 of meer"] <- 9
table(all_data_forstata$income)


all_data_forstata$interest <- NA
all_data_forstata$interest[all_data_forstata$QPT3.1 == "Helemaal niet geÃ¯nteresseerd"] <- 0
all_data_forstata$interest[all_data_forstata$QPT3.1 == "niet geÃ¯nteresseerd"] <- 1
all_data_forstata$interest[all_data_forstata$QPT3.1 == "Nauwelijks geÃ¯nteresseerd"] <- 2
all_data_forstata$interest[all_data_forstata$QPT3.1 == "Tamelijk geÃ¯nteresseerd"] <- 3
all_data_forstata$interest[all_data_forstata$QPT3.1 == "Heel erg geÃ¯nteresseerd"] <- 4
table(all_data_forstata$interest)

all_data_forstata$leftright <- NA
all_data_forstata$leftright[all_data_forstata$QPT3.4 == "Weet ik niet/wil ik niet zeggen"] <- NA
all_data_forstata$leftright[all_data_forstata$QPT3.4 == "0 Links"] <- 0
all_data_forstata$leftright[all_data_forstata$QPT3.4 == "1"] <- 1
all_data_forstata$leftright[all_data_forstata$QPT3.4 == "2"] <- 2
all_data_forstata$leftright[all_data_forstata$QPT3.4 == "3"] <- 3
all_data_forstata$leftright[all_data_forstata$QPT3.4 == "4"] <- 4
all_data_forstata$leftright[all_data_forstata$QPT3.4 == "5"] <- 5
all_data_forstata$leftright[all_data_forstata$QPT3.4 == "6"] <- 6
all_data_forstata$leftright[all_data_forstata$QPT3.4 == "7"] <- 7
all_data_forstata$leftright[all_data_forstata$QPT3.4 == "8"] <- 8
all_data_forstata$leftright[all_data_forstata$QPT3.4 == "9"] <- 9
all_data_forstata$leftright[all_data_forstata$QPT3.4 == "10 Rechts"] <- 10 
table(all_data_forstata$leftright)

all_data_forstata$risks <- NA
all_data_forstata$risks[all_data_forstata$QPT3.5 == "Weet ik niet"] <- NA
all_data_forstata$risks[all_data_forstata$QPT3.5 == "Helemaal niet"] <- 0
all_data_forstata$risks[all_data_forstata$QPT3.5 == "Niet"] <- 1
all_data_forstata$risks[all_data_forstata$QPT3.5 == "Neutraal"] <- 2
all_data_forstata$risks[all_data_forstata$QPT3.5 == "Graag"] <- 3
all_data_forstata$risks[all_data_forstata$QPT3.5 == "Heel graag"] <- 4
table(all_data_forstata$risks)

all_data_forstata$salience_char <- NA
all_data_forstata$salience_char[all_data_forstata$salience == 0] <- "zero"
all_data_forstata$salience_char[all_data_forstata$salience == 1] <- "one"
all_data_forstata$salience_char[all_data_forstata$salience == 2] <- "two"
all_data_forstata$salience_char[all_data_forstata$salience == 3] <- "three"
all_data_forstata$salience_char[all_data_forstata$salience == 4] <- "four"
all_data_forstata$salience_char[all_data_forstata$salience == 5] <- "four"
all_data_forstata$salience_char[all_data_forstata$salience == 6] <- "four"
table(all_data_forstata$salience_char)

all_data_forstata$salience_char <- as.factor(all_data_forstata$salience_char)
all_data_forstata$salience_char <- relevel(all_data_forstata$salience_char, "one")







all_data_stata <- all_data_forstata %>% dplyr::select(rid, salience, QPT2.1 , QPT2.2 , QPT2.3 , QPT2.4 , QPT3.1 , 
                                      QPT3.2 , QPT3.4 , QPT3.5, age, gender, education, income,
                                      interest, leftright, risks, salience_char, issue)
library(foreign)
write.dta(all_data_stata, "data/all_data_long.dta") # Save the data for STATA, the stata file then produces the coefficients, which have then been formatted in LaTeX









##### Check distribution of attributes (Figure A2) #####
library(cregg)

dd <- all_data
dd$strategy_recoded <- as.factor(dd$strategy_recoded)
dd$Partij <- as.factor(dd$Partij)
dd$issue <- as.factor(dd$issue)

f1 <- selected~strategy_recoded + Partij + issue
plot(cj_freqs(dd, f1, id = ~rid))












###### Number of votes by issue
#number of votes by issue #Figure 5, right panel
#all_data_fordes <- readRDS("data/final_data_wide_all_recoded.rds")
#all_data_fordes <- all_data_fordes %>% dplyr::select(rid, QV_eu, QV_climate, therm_pvda,
#                                                     QV_farms, QV_lgbt, QV_immi, QV_wage) %>% 
#  distinct(rid, .keep_all=T) %>%
#  pivot_longer(cols = starts_with("QV_"),
#               names_to="Issue", names_prefix = "QV_",
#               values_to = "Votes") %>% dplyr::select(-rid)
#
# ##Save only the version with exact the variables used (JoP reproduction requirement)
#write_rds(file="data/final_data_wide_all_recoded_reprod.rds", x=all_data_fordes)
all_data_fordes <- readRDS("data/final_data_wide_all_recoded_reprod.rds")


all_data_fordes$Issue[all_data_fordes$Issue=="eu"] <- "EU"
all_data_fordes$Issue[all_data_fordes$Issue=="climate"] <- "Climate"
all_data_fordes$Issue[all_data_fordes$Issue=="farms"] <- "Farms"
all_data_fordes$Issue[all_data_fordes$Issue=="lgbt"] <- "LGBTQ+"
all_data_fordes$Issue[all_data_fordes$Issue=="immi"] <- "Immigration"
all_data_fordes$Issue[all_data_fordes$Issue=="wage"] <- "Min. Wage"
all_data_fordes$Issue <- as.factor(all_data_fordes$Issue)
levels <- c("Climate", "EU", "Farms", "Immigration", "LGBTQ+", "Min. Wage")
all_data_fordes$Issue <- factor(all_data_fordes$Issue, levels = levels)
all_data_fordes$Issue <- factor(all_data_fordes$Issue, levels = rev(levels(all_data_fordes$Issue)))

library(ggridges)
cols <- c("#760B08", "#B51811", "white", "#469C40", "#347731")
values <- c(0,0.2,0.5,0.7, 1)

votes_overview <- ggplot(data=all_data_fordes, aes(x=Votes, y=Issue, fill=stat(x))) +  
  geom_density_ridges_gradient(size=0.8, bandwidth=0.37, scale = 0.94) + scale_fill_gradientn(colors = cols, values = values) + 
  labs(y="Density", x="Votes") +  ggtitle("", subtitle="Distribution of QV votes") +
  theme(axis.text.y=element_blank(), legend.position = "none") +
  scale_x_continuous(limits = c(-6, 6), breaks = c(-6:6), 
                     labels = c("-6 Anti", "-5", "-4", "-3", "-2", "-1", 
                                "0", "1", "2", "3", "4", "5", "6 Pro")) +

  scale_y_discrete(expand = expansion(add = c(0.1, 0.8)),  position = "right")


votes_overview #Figure 5, right panel.  963
#ggsave("plots/fig5r_votes_overview2.png", plot=votes_overview, dpi=600, width=7.29, height = 7.05)
#rm(cols, values, all_data_fordes)

fig5both <- ggarrange(simulation_counterfactuals, votes_overview, nrow=1)
ggsave("plots/fig5both.png", plot=fig5both, dpi=600, width=11, height = 7)







#### Wide plot for partisanship without issue salience (Figure A3, left panel)
# These are models without issue salience
just_p_elec <- (lm_robust(selected~strategy_recoded * partisan_altelect + issue + 
                            Partij * strategy_recoded, data=all_data, clusters = rid))
#only_p_margins_elec <- margins_summary(just_p_elec, at=list(partisan_altelect=c(0,1)))

just_p_partisan <- (lm_robust(selected~strategy_recoded * partisan + issue + 
                                Partij * strategy_recoded, data=all_data, clusters = rid))
#only_p_margins_partisan <- margins_summary(just_p_partisan, at=list(partisan=c(0,1)))
#only_p_margins_partisan <- only_p_margins_partisan %>% rename(Strategy = factor, Partisanship=partisan) %>%
#  mutate(Estimate="Election today & Thermostat")

partisanonly <- only_p_margins_elec
partisanonly <- partisanonly %>% rename(Strategy = factor, Partisanship = partisan_altelect) %>%
  mutate(Estimate="Election today")

partisanonly <- full_join(partisanonly, only_p_margins_partisan )
partisanonly$Partisanship <- as.factor(partisanonly$Partisanship)

#Create the plot
partisan_only <- ggplot(data=subset(partisanonly, Strategy == "strategy_recodedAmbiguity"), 
                        aes(x=Partisanship, y=AME, ymin=lower, ymax=upper, color=Estimate, shape=Estimate)) +
  geom_pointrange(position=position_dodge(width=0.2)) + 
  labs(x="", y="AMCE", subtitle="The Effect of Ambiguity as Opposed to Agreeing by Partisanship") +
  scale_colour_grey(end = 0.6) +
  scale_x_discrete(labels=c("Non-Partisan", "Partisan")) + theme(legend.position = "bottom")
partisan_only # Figure A3, left panel
#ggsave("plots/figa3L_partisan_only.png", plot= partisan_only, dpi=600, width=3584, height=2616, units = c("px"))

figa3both <- ggarrange(partisan_only, with_p_margins_plot_both, nrow=1)
figa3both
ggsave("plots/figa3both.png", plot=figa3both, dpi=600, width=12, height = 7)
ggsave("plots/figa3both.tiff", plot=figa3both, dpi=600, width=12, height = 7, device="tiff")





# And save both panels for figure 1
fig1both <- ggarrange(without_p_margins_plot_both, without_p_predp_plot, nrow=1)
ggsave("plots/fig1both.png", plot=fig1both, dpi=600, width=7168, height=2616, units = c("px"))
ggsave("plots/fig1both.tiff", plot=fig1both, dpi=600, device="tiff", width=7168, height=2616, units = c("px"))





## Save the objects that take longer to run so that it is easier to replicate the figures
library(gdata)
gdata::keep(only_p_margins_partisan, only_p_margins_elec, without_p_margins_pos, with_p_margins_0, with_p_margins_1, 
     without_p_margins, only_p_margins_partisan, only_p_margins_elec, sure=T)
